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ABSTRACT 

We study accretion processes for tidally disrupted stars approaching supermassive 
black holes on bound orbits, by performing three dimensional Smoothed Particle Hy- 
drodynamics simulations with a pseudo-Newtonian potential. We find that there is 
a critical value of the orbital eccentricity below which all the stellar debris remains 
bound to the black hole. For high but sub-critical eccentricities, all the stellar mass 
is accreted onto the black hole in a finite time, causing a significant deviation from 
the canonical £ -5 / 3 mass fallback rate. When a star is on a moderately eccentric orbit 
and its pericenter distance is deeply inside the tidal disruption radius, there can be 
several orbit crossings of the debris streams due to relativistic precession. This dissi- 
pates orbital energy in shocks, allowing for rapid circularization of the debris streams 
and formation of an accretion disk. The resultant accretion rate greatly exceeds the 
Eddington rate and differs strongly from the canonical rate of t -5 / 3 . By contrast, there 
is little dissipation due to orbital crossings for the equivalent simulation with a purely 
Newtonian potential. This shows that general relativistic precession is crucial for ac- 
cretion disk formation via circularization of stellar debris from stars on moderately 
eccentric orbits. 

Key words: accretion, accretion discs - black hole physics - gravitational waves - 
galactic: nuclei - hydrodynamics 
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1 INTRODUCTION 

Supermassive black holes (SMBHs) larger than 10 5 Mq are 
now known to exist ubiquitously in galactic nuclei. SMBHs 
in nearby galaxies can be studied dynamically, whereas in 
more distant galaxies only the ~ 1% of SMBHs undergoing 
major accretion episodes can be easily observed. The tidal 
disruption and subsequent accretion of a star by a SMBH, 
although intrinsically a rare event, is of observational in- 
terest as a way to probe dormant SMBHs in the local uni- 
verse, because it can produ ce a powerful flare in excess of 



a seed black hole grows into a SMBH (|Zhao et aL []2002l : 



the Eddington luminosity (iCarter & Luminet _ 19831; Reed 



gity ^ 

Il988l : lEvans fc Kochanekll 19891 : IStrubbe fc Quataertll2009L 
These events are also of interest to low-frequency gravita- 
tional wave astronomy, as the stellar dynamical processes 
which funnel stars into the low angular momentum orbits 
necessary for tidal disruption events (TDEs) are similar to 
those which produce the extreme-mass rati o inspiral of a 
stella r -mass compact object onto a SMBH ([Frank &: Reed 
Il976l : IWang fe Merrittl 120041 : iMadigan et all l201lh . TDEs 
have also been considered as part of the means by which 
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iMiralda-Escude fc Kollmeierll2005l : lBromlev et alJl2012h . All 

of these aspects of stellar tidal disruption have motivated 
past studies of TDEs. 

The traditional picture of a TDE proceeds as follows: 
a star at large separation (~ 1 pc) approaches a SMBH on 
a nearly parabolic orbit. After the star is tidally disrupted 
by the black hole, half the stellar debris becomes gravita- 
tionally bound to the SMBH, because it loses orbital energy 
inside the tidal radius. The bound debris falls back, and, 
after circularizing due to collisional shocks with other gas 
streams, accretes onto the black hole. Kepler's third law im- 
plies that the mass return rate d ecays with a -5/3 power of 
time (lReedll988l : [Phinnevlll989h . asymptotically approach- 
ing zero. This simple analytic picture has been validated to 
the fi rst order of approximation by hydrodynamical simula- 
tions ([Evans fc Kochan ek 1989) , albeit with deviations from 
this law at early times ([Lodato et al.ll2009h . Similar power- 
law behavior for the flare lightcurve is often assumed (i.e. 
L oc M\ although her e the theoretical evidence is less clear 
(lLodato fc Rossil201lh . 

All-sky surveys in the X-ray and UV have so 
far observed 13 candidate tidal disruption flares 
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Komossa fe Badd 119991: iGruue e t al. 19991: Ig reiner et al.l 
2000t [ Komossa et al II2QQ4J: iHalpern et al.|[20o3:lDogiel et al. 1 
2QQ9: lMaksvm et alfeoToUBloom et al.ll201ll : lBurrows et al.1 



39; 
11) 



20 111 ). The observed light curves are in reasonable agree- 



ment with the theoretically predicted mass fallback rate o f 



t 



-5/3 



although some show deviations ([Burrows et al.|[201ll ) 
and the number of samples is sufficiently small to make 
detailed testing of theoretical models difficult. A smaller 
number of stron g TDE candidates have been found at opt i- 
cal wavelengths flvan Velzen et al ][20ll iGezari et alJ lioi^ . 
Notably, two of the best-sampled TDEs differ strongly from 
the canonical theoretical picture: one posse sses relativistic 
jets ([Bloom et al.ll201ll : IZauderer et al.ll201lh. and the other 
lacks hydrogen lines in its spectra (|Gezari et al.l l2012h . 

It has been inferred from observations that the 
event rate of tidal disruption is ~ 10 _5 yr _1 per galaxy 
([Donley et al.l l2002h . This observed rate is in rough agree- 
ment with uncertain theoretical rate estimates. In the stan- 
dard scenario, stars are supplied to the SMBH by two body 
scattering, on a stellar relaxational timescale. Stars with an- 
gular momentum less than a critical value are in the phase 
space "loss cone" and are tidally destroyed on a dynami- 
cal time. Theoretical calculations indicate this rate to be 
~ 10" 4 - 10~ 6 yr" 1 for Milky Way-like galaxies, and that 
the peak flux into the loss cone comes from parsec scales. 
This motivates the assumption of nearly parabolic orbits 
([Magorrian fc Tremainelll999l : IWang fc Merrittll2004h . 

However, there are alternate sources of TDEs distinct 
from the standard two-body scattering model. Many of these 
feed stars to the SMBH at lower eccentricity. Our aim here is 
to quantify through hydrodynamical numerical simulations 
how the observable properties of tidal disruptions of stars 
on eccentric orbits deviate from the standard ones. In § [2] 
we describe numerical procedures for simulating TDEs. We 
then analyze results of our simulations in § [3] In § [4] we 
consider nonstandard sources of TDEs and whether they can 
supply stars on low eccentricity orbits to a SMBH. Finally, 
§ [5] is devoted to summary and discussion of our scenario. 



2 NUMERICAL METHOD 

In this section, we describe our procedures for numerically 
modeling the tidal disruption of stars on bound orbits. The 
simulations presented below were performed with a three- 
dimensional (3D) SPH code, which is a particle method that 
divides the fluid into a set of discrete "fluid elements" (i.e. 
particles) , and is flexible in setting various initial configura- 
tions. Th£Code_jsJ^ by 
iBenz I (ll990h : lBenz et~al~l (|l990T ): lBate et al.1 dl995h and has 
been e xtensively used by ma ny papers (e.g., lOkazaki et al.l 
(|2002MHavasaki et al.1 ([20071 )1. The SPH equations with the 
standard cubic-spline kernel are integrated using a second- 
order Runge-Kutta-Fehl berg integrator w ith individual time 
steps for each particle ([Bate et al.l IT995), which results in 
enormous computational time savings when a large range of 
dynamical timescales are involved. In simulations shown in 
this paper, we adopt standard SPH artificial viscosity pa- 
rameters aspH = 1 or /3sph = 2, unless otherwise noted. We 
have performed two-stage simulations. We model a star as a 
polytropic gas sphere in hydrostatic equilibrium. The tidal 



disruption process is modeled by setting the star in motion 
through the gravitational field of an SMBH. 

2.1 Formation of polytropes 

In our code, the polytrope is composed of an ensemble of gas 
particles, each of which has a mass chosen to be 10 -5 M© 
with a variable smoothing length. The particles are initially 
uniformly distributed in a spherical fashion, with an ini- 
tial temperature Ti n i = 1.2 X 10 6 i^. The initial spherical 
gas cloud is allowed to collapse under self-gravity using the 
polytropic equation of state: 

P = Kp 1+1/n , (1) 

where n is the polytropic index and K is assumed to be 
kept constant throughout the collapse. Note that n = 1.5 
corresponds to 7 = 5/3. The simulations continue over five 
dynamical times, where the dynamical time is defined by 

•(£)""(£)■%. m 

where G is Newton's constant and r* and ra* are the stellar 
radius and mass, respectively. 

Fig. [1] shows the radial density profile of the polytropic 
gas sphere at 4.4Q -1 , where the magnitude of ratio of ther- 
mal and gravitational energy becomes ~ 0.5. While circle 
marks show the radial density profile obtained from the SPH 
simulations, the solid line shows the profile obtained by nu- 
merically solving the Lane-Emden equation: 



(3) 



where we define as dimensionless quantities, £ = r/r c with 

and n = p/pc with the central density p c . With boundary 
conditions that — 1 and dO/dt; = at £ = and — and 
£ = Rq/tc at the surface of the star for 1 solar mass sun- 
type star, we obtained that p c ~ 8.4 g cm -3 , r c — 0.27 Rq 
and K-2.5x 10 15 . 



2.2 Initial configuration 

A star is tidally disrupted when the tidal force of the black 
hole acting on the star is stronger than the self- gravitational 
force of the star. The radius where these two forces balance 
is defined as the tidal disruption radius by 



/M BH \ 1/3 
n r*, 



(5) 



where Mbh is the black hole mass. 

Fig. [2] shows a density map of the polytropic gas sphere 
at 4.4Q" 1 . The left panel shows the column density of poly- 
tropic sphere with n = 1.5 over a range of five orders of mag- 
nitude. The right panel shows the star-black hole system on 
the x-y plane, where both axes are normalized by rt and 
the black hole is put on the origin of the system. The initial 
position of the star is given by ro = (ro cos 0o, ro sin 0o, 0), 
where |ro| = To is the radial distance from the black hole 
and 0o shows the angle between x-axis and ro- The star 
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Figure 1. The radial density distribution of the polytropic gas 
sphere with n = 1.5. The solid line and circle marks show the 
solution for the Lane-Emden equation with n = 1.5 and the solu- 
tion derived from the SPH simulation, respectively. The density is 
normalized by the central density p c = 8.4 gem -3 and £ = r/r c , 
where r c ~ 0.27 Rq. 



within a small square is zoomed out to the whole left panel. 
Note that the stellar size on the plot has been scaled up by 
a factor of a few for visual clarity. 

The motion of a test particle in the central SMBH po- 
tential U(r) admits two conserved quantities. The radial ve- 
locity and angular velocity are given by e nergy conservation 
and a ngular momentum conservation as ([Landau fc Lifshitzl 
Il969h 



(6) 
(7) 



where e and I are the specific energy and the specific angular 
momentum, respectively. For bound orbits, 



(r p /r a ) 2 t/(r p ) - t/(r a ) 
(rpA-a) 2 - 1 



I = y/2rl(e ~ U(r p )) = ^2ri(e - £7(r a )), 



(8) 
(9) 



where r p and r a are the pericenter distance and 
the apocenter distance, respectively. Therefore, the ini- 
tial velocity vector is given by Vo = (r(ro) cos 0o — 
r o 0(0o) sin O , r(r ) sin O + r o 0(0o) cos 0o, 0). 

In our simulations, the black hole is represented by a 
sink particle with the appropriate gravitational mass Mbh- 
All gas particles that fall within a specified accretion radius 
are accreted by the sink particle. We set the accretion radius 
of the black hole as equal to the Schwarzshild radius rs = 
2GMbh/c 2 , with c being the speed of light. 



2.3 Treatment of relativistic precession 

In order to treat approximately the relativistic precession of 
a test particle in the Schwarzschild metric, we incorporate 
into our SPH code the following pseudo-Newtonian potential 
(|Wegdl2012h : 



U(r) 



GM BH 



ci + 



1 



Cl 



l-c 2 (Vs/2r) 2r 



(10) 



(-4/3)(2 + V6), c 2 = (4^6-9), 



where we adopt that c\ 
and c 3 = (-4/3)(2V6 - 3). Equation JTJJ reduces to the 
Newtonian potential when c\ — 1 and c 2 = C3 =0, while 
it reduces to the wel l-known Paczynski fc W iita pseudo- 
Newtonian potential ([Paczviiski fc Wiitalll980h when ci = 
C3 = 0. Note that equation (|10p includes no higher-order 
relativistic effects such as the black hole spin or gravitational 
wave emission. 

Fig. [3] tests how accurately the SPH particles compos- 
ing the star follow test particle motion in the Schwarzschild 
metric. Specifically, the figure shows the motion of three 
particles orbiting the black hole with (e, /3) = (0.8, 5), where 
f3 = n/r p is the penetration factor which determines how 
deeply the star plunges into the black hole potential inside 
the tidal disruption radius. The central point and dashed 
circle show the black hole and the tidal disruption radius, 
respectively. The small crosses, dotted line, and solid line 
represent the geodesies of the Schwarzschild metric, orbits 
of the test particle under the pseudo-Newtonian potential, 
and those of a SPH particle in the simulation of Model 2a, 
respectively. 

During the first one and a half orbits, the test parti- 
cle is in good agreement with the corresponding geodesic. 
The subsequent deviation from the geodesic is due to the 
relatively low eccentricity and high f3, sinc e equation (JTO)) is 
tailored to parabolic orbits with lower (3 ([Wegdl2012n (see 
also Fig. 2}. The first three and half orbits of the test par- 
ticle are in rough agreement with corresponding orbits of 
the SPH particle. Subsequent large deviations originate from 
orbital circularization via shock-induced energy dissipation 
(see Section \3. 2 1 in detail), and although our model becomes 
unreliable at this point, it is clear that stream crossing is 
leading to rapid energy dissipation. 

Fig. [4] shows the dependence of the error rate on the 
eccentricity of a test particle moving under the pseudo- 
Newtonian potential given by equation (|10|) . The error rate 
is defined by (A0gr — A0pn)/A0gr x 100, where A0gr 
and A0pn are the precess i on an gles corresponding to equa- 
tion (2) and (4) of lWegd (|2012h . respectively, with a given 
specific angular momentum, specific orbital energy, and 
equation ([TO]) . The solid line, dashed line, and dotted line 
show error rates of an orbit with r p = 3rs (/3 ~ 7.9), 
r p = 5rs (/3 ~ 4.7), and r p = 10rs (fi ~ 2.4), respectively, 
where r t ~ 23.6rs is adopted. The error rate increases as the 
eccentricity is lower and the penetration factor is higher. 
From the figure, the error rate is estimated to be roughly 
15% for e = 0.8 with f3 = 5, while it is less than 1% for 
e > 0.95 with arbitrary /3. 



2.4 Numerical models 

We have performed eight simulations of tidal disruption 
events with different parameters. The common parame- 
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Figure 2. Initial configuration of our simulations. The left panel shows a density map of the polytropic gas sphere with n = 1.5. Its radius 
is IRq and the number of SPH particles iVspH are annotated in the right-bottom corner of the panel. The panel on the right-hand-side 
shows the initial position of the star on the x-y plane. The small square corresponds to the entire left panel, the dashed circle shows the 
tidal disruption radius r*t and the black hole is positioned at the center. The simulation run time t in the right-top corner is normalized 
by Q* 1 . In the right panel, the star is initially located at (ro cos <po, ro sin^o), where ro = 3rt and <po = — OAtt. 



Table 1. The first column shows each simulated scenario. The second, third, fourth, and fifth columns are the penetration factor 
(3 = r p /rx, the initial orbital eccentricity e*, the initial semi-major axis a*, and the radial distance between the black hole and the initial 
position of the star, respectively. The last column describes the remark for each model. 



Model f3 = rt/r p e* a* [n] ro [n] Remarks 



la 


1 


1.0 




3 


Newtonian 


lb 


1 


1.0 




3 


Pseudo-Newtonian 


lc 


5 


1.0 




3 


Newtonian 


Id 


5 


1.0 




3 


Pseudo-Newtonian 


2a 


5 


0.8 


1.0 


1.8 


Pseudo-Newtonian 


2b 


5 


0.8 


1.0 


1.8 


Newtonian 


3a 


1 


0.98 


50.0 


3 


Pseudo-Newtonian 


3b 


5 


0.98 


10.0 


3 


Pseudo-Newtonian 



ters through all of simulations are following: ra* = 1M©, 
r* = IRq, M BH = 1O 6 M , O = -OAtt, and 7 = 5/3. The 
total number of SPH particles used in each simulation are 
10 5 , and the termination time of each simulation is 
Table 1 summarizes each model. Model la shows the stan- 
dard TDE under the Newtonian potential, while Model lb 
has the same simulation parameters as Model la, except 
that the star moves under the pseudo-Newtonian potential 
given by equation ([TO]) . Model lc and Id have the same pa- 
rameters as Model la and Modellb, respectively, but for 
/3 = 5. Model 2a has the same parameters as Model Id ex- 
cept that the star is on an eccentric orbit, with e = 0.8. 
Model 2b is the same parameters as Model 2a except that 
the star moves under the Newtonian potential. Model 3a 
has the same parameters as Model 2a but for e = 0.98 and 
/3 = 1, and Model 3b has the same parameters as Model 3a 
but for = 5. 



3 TIDAL DISRUPTION OF STARS ON 
BOUND ORBITS 

We first describe the evolution of a tidally disrupted star 
for standard TDEs (e = 1). Next, the process of tidal dis- 
ruption for a star on a fairly eccentric orbit (e = 0.98) is 
presented. Finally, accretion disk formation in relatively low 
eccentricity TDEs (e = 0.8) is presented. 

As an approaching star enters into the tidal disruption 
radius, its fluid elements become dominated by the tidal 
force of the black hole, while their own self- gravity and pres- 
sure forces become relatively negligible. The tidal force then 
produces a spread in specific energy of the stellar debris 



The total mass of the stellar debris is denned with the dif- 
ferential mass distribution m(e) = dM(e)/de by 
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Figure 3. Motion of three particles with orbital parameters 
(e, (3) = (0.8,5) in the x-y plane. Each axis is normalized by 
r*t given by equation (|5)>. The central point and dashed circle 
represent the black hole and tidal disruption radius, respectively. 
The small crosses, dotted line, and solid line show the geodesies 
in the Schwarzshild metric, orbits of a test particle under the 
pseudo-Newtonian potential given by equation (jlpp . and orbits of 
a SPH particle in the simulation of Model 2a (see text of § 13. 2$ , 
respectively. 
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Figure 4. Dependence of error rate of precession angle A<£ on the 
eccentricity of a test particle moving under the pseudo-Newtonian 
potential given by equation (|10|) . The solid line, dashed line, and 



dotted line show error rates of an orbit with r D 



: 3r s , r p = 5r s , 



and r p = 10rs, respectively. The horizontal dash-dotted line and 
dash-three-dotted line show the error cutoffs of 10% and 20%, 
respectively. 



M(e) 



J — c 



i(e )de 



(12) 



When a star is disrupted from a parabolic orbit, m(e) will 
be centered on zero and distributed over — Ae ^ e ^ Ae. 

Since the stellar debris with negative specific energy 
is bounded to the black hole, it returns to pericenter and 
will eventually accrete onto the black hole. If its specific 
energy is approximately equal to the specific binding en- 
ergy, e ~ — GMbh/^cl with the semi-major axis of the stel- 
lar debris a, the mass fall back rate is then defined by (e.g., 
lEvans fe Kochanekl (| 1989ft ) 



dM(e) 



de 



(e<0), 



dM 
dt 
where 
de _ 1 
dt ~3 ( 

This is derived from the relationship between the orbital 
period and the specific energy: 

7T GMbh 



-(27TGMBH 



2/3^-5/3 



(13) 



(14) 



V2 e 3/2 



(15) 



from Kepler's third law. The standard dM /dt oc £ _5//3 scal- 
ing then arises with the assumption that dM/de is constant, 
or at least asymptotes to a constant value at late times. 
From equation (|15[) , the orbital period of the most tightly 
bound orbit (t m in) and that of the most loosely bound orbit 
(tmax) are given by 



7T 

/a/2) (Mbh/ 

= OO 



\l/2 



(16) 



(e = Ae) 
(e = 0) 

Figs. [5]|8]show the evolution of differential mass distri- 
butions and their corresponding mass fallback rates in Mod- 
els la- Id. While the differential mass distribution is shown 
in panel (a), the mass fallback rate is shown in panel (b). 
In both panels, the solid line, dashed line, dot-dashed line, 
and dotted line represent the mass distributions and corre- 
sponding mass fallback rates at t = 0.0, t = 0.37, t = 0.79, 
and t = 4.0, respectively. Note that t = 0.5 corresponds to 
the time when the star reaches pericenter. In panel (b), the 
horizontal solid line denotes the Eddington rate: 

1 ^Edd 

rj c 2 

M BH 



2.2 x 10" 



(£/'( 



V 10 6 M G 



M yr _ 



(17) 



where Lsdd = 47vGMBB\m p c/<JT is the Eddington luminos- 
ity with m p and <jt denoting the proton mass and Thomson 
scattering cross section, respectively, and 77 is the mass-to- 
energy conversion efficiency, which is set to 0.1 in the fol- 
lowing discussion. 

In Model la, the mass distribution broadens with time 
as the stellar debris is put on near-ballistic orbits by tidal 
interaction with the black hole. In panel (a), the central 
peak formed after t > 1.0 is due to mass congregation, from 
the self-gravity of the stellar debris (in this scenario, the 
star is barely disrupted). The energy spread corresponds to 
Ae before and after the tidal disruption. The correspond- 
ing mass fallback rates are proportional to t _5//3 except for 
the solid line in panel (b) . This is in good agreemen t with 
the literatures dReesI Il988l : lEvans fc Kochanekl [l989h . The 
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slight deviation from t to the —5/3 power originates from 
the convexity around Ae d ue to re-congregatio n of stellar 
mass under self- gravity fsee lLodato et al l ([2009h l. 

Model lb has the same simulation parameters as 
Model la except that it adopts the pseudo-Newtonian po- 
tential. In panel (a) of Fig. [6] the final energy spread and 
the central peak are slightly wider and milder than those of 
Model la. This is because the energy imparted by the tidal 
force under the pseudo-Newtonian potential is slightly larger 
than that of the Newtonian potential. From equation (|10[) , 
the energy spread under the pseudo- Newtonian potential is 
evaluated to be 



AepN 



dU(r) 



dr 



r* = Ae 



1 — ci rs 
Cl + [l-c 2 (rs/2r t )P + 2 ° 3 2^ 



The expected energy spread AepN ~ 1.1 Ae for 1O 6 M0 is in 
good agreement with the simulated energy spread. 

Model lc has the same simulation parameters as 
Model la except for (3 = 5, while Model Id has the 
same simulation parameters as Model lc except for adopt- 
ing the pseudo-Newtonian potential. For higher value of 
P than unity, the tidal forces acting on the star become 
stronger because the pericenter distance is smaller than 
the tidal disruption radiu s. According to equation (2) of 
lEvans fe Kochanekl (| 19891 ). the energy spread of the stel- 
lar debris can be estimated to be /3 2 Ae, and should become 
25Ae for Models lc and Id. In these two models, however, 
the energy spread is only about 1.4 times wider than that 
of standard model, showing that the energy spread weakly 
depe nds on (3. This is in agreement with recent numer- 
ical^ Guillochon Ramirez- Ruizl l2012h and analytic work 
(jStone et al. 1 12012) that demonstrates the inaccuracy of the 
traditional formula for Ae, in place of which we have used 
the more accurate equation (|11[) . Since the tidal disruption is 
less marginal as /3 is higher, the re-congregation of the mass 
due to self- gravity of the stellar debris is prevented. This 
leads to the mildly-sloped mass distribution, and therefore 
the peak of the mass fallback rate also smooths. There is 
no remarkable difference between Model lc and Model Id 
except that the energy spread of Model Id is slightly wider 
than that of Model lc following equation JTJ 



3.1 Eccentric tidal disruption 

The specific orbital energy of a star on an eccentric orbit is 
given by 



Corb 



GMbk 
2a* 



GM B 
2r t 



~P(1 



(19) 



This quantity is less than zero because of the finite value of 
a*, in contrast to the standard, parabolic orbit of a star. If 
e or b is less than Ae, all the stellar debris should be bounded 
to the black hole, even after the tidal disruption. The con- 
dition e or b = Ae therefore gives a critical value of orbital 
eccentricity of the star 

1/3 

ecrit « 1 - - ( ^P- ) , (20) 



2 / m* y 



below which all the stellar debri s should remain gravitation - 
ally bound to the black hole (|Amaro-Seoane et al.l 1201 lh . 
The critical eccentricity is evaluated to be e cr it = 0.996 for 
Model 2a and Model 3b, whereas e cr it = 0.98 for Model 3a. 



For the eccentric TDEs, t m i n and t max are obtained by 
substituting e = Ae ± e or b into equation JT5J). The orbital 
period of the most tightly bound orbit is given by 



nz 1 



s/2P 3/2 {l -e«) 3 / 2 ' 



(21) 



The orbital period of the most loosely bound orbit t-max IS 
estimated to be oo for e* ^ e cr it, while t max converges with 



V2 



n* 1 



P(l 



1/3 



1 -3/2 



VMb 



(22) 



for e* < e cr it- The duration time of mass fallback for eccen- 
ic TDEs with e* < e cr it is thus predicted to be finite and 




(23) 



Evaluating this gives At « 4.3Q" 1 for Model 2a and At « 
207Q" 1 for Model 3b, whereas At ^ oo for Model 3a. 

Figs. [9] and [10] show the evolution of differential mass 
distributions and corresponding mass fallback rates in 
Model 3a and Model 3b, respectively. The figure formats are 
the same as Figs. [5][8] The mass is not distributed around 
zero but around — Ae in Model 3a, and around — 5Ae in 
Model 3b. This is because the specific energy of the initial 
stellar orbit is originally negative (see equation I19[) . It is 
clear from the negative shift of the mass distribution's cen- 
ter that most of the mass in Model 3a and all of the mass 
in Model 3b are bounded. As shown in panel (a) of Fig. [9] 
the resultant energy spread is slightly larger than we analyt- 
ically expected. This suggests that the critical eccentricity 
is slightly smaller than the value in equation (|20p . 

From equations (|21[) and (|22[) , t m in = 0.04yr and t max —> 
oo for Model 3a, and t m i n = 0.0035yr and t max = 0.014yr 
for Model 3b. These values deviate from simulation results, 
because they were derived assuming a smaller spread in en- 
ergy than in our simulations. Furthermore, t m i n is shorter 
than that of the parabolic case. This makes the mass fall- 
back rate about one order of magnitude higher. Notably, the 
mass accretion is completely finite in Model 3b and its rate 
is enhanced by about two orders of magnitude, because its 
timescale is much shorter than that of the standard model. 



3.2 Accretion disk formation 

Here we describe how an accretion disk forms rapidly around 
the black hole in Model 2a. The orbital angular momentum 
of a star passing inside the tidal disruption radius should be 
conserved before and after the tidal disruption, if there is 
no mechanism to redistribute the angular momentum. The 
orbits of the stellar debris can then be circularized by dissi- 
pation of orbital energy, primarily due to shocks from orbit 
crossing, which conserve orbital angular momentum. Debris 
semi-major axes will thus approach the circularization ra- 
dius of the initial stellar orbit. 

The specific energy and specific angular momentum of 
a test particle with the orbital parameters of Model 2a, 
moving under the pseudo-Newtonian potential, are given as 
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Figure 5. Evolution of differential mass distributions over specific energy of stellar debris and their corresponding mass-fall back rates, 
in Model la. In panel (a), the dashed line (red), dash-dotted (fairly green), dotted line (blue) and solid line (black) show the differential 
mass distributions measured at t = 0, t = 0.37, t = 0.79, and t = 4, respectively. The energy is measured in units of Ae given by 
equation JTTJ. In panel (b), the dash line (red), dashed-dotted line (fairly green), dotted line (blue), and solid line (black) show the 
mass fall-back rates measured at t = 0, t = 0.37, t = 0.79, and t = 4, respectively. Each mass fallback rate is evaluated by using each 
differential mass distribution represented in panel (a). The mass fallback rate and the orbital period are normalized by 1 solar mass per 
year and a year, respectively. 
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Figure 8. Same as Fig. [5] but for Model Id. 




et P = —0.497 and h P = 0.714, respectively, by substituting 
(r p ,r a ) = (0.2rt, 1.8rt) with (e,/3) = (0.8,5) into equations 
(|9]) and dTOj) - Fig. HH shows the evolution of the specific en- 
ergy, specific internal energy, and the specific angular mo- 
mentum in Model 2a. In panel (a), the solid line and the 
dashed line show the SPH specific energy and that of a test 
particle, respectively. The small inside panel shows the evo- 
lution of the specific internal energy, which is dissipated and 
radiated away after t = 2.4Q -1 due to the shocks from self- 
intersections of the debris orbits induced by relativistic pre- 
cession. Since the internal energy contributes negligibly to 
the total energy budget, the specific energy approximately 
equals the specific binding energy. On the other hand, panel 
(b) shows the evolution of the specific angular momentum 
normalized by lt p . Its marginal fluctuation is a numerical er- 
ror with a magnitude of 0.15%, indicating that the specific 
angular momentum is almost conserved. Assuming that the 
circularized disk has an axisymmetric Keplerian rotation, we 
can evaluate that the circularization radius r c irc is ~ 2.5r p . 

Fig.[T2lshows sequential snapshots of the surface density 
of stellar debris (projected on the x-y plane in five orders of 
magnitude, in a logarithmic scale) for Model 2a. The central 
small point, dashed circle, and dotted line show the black 
hole, tidal disruption radius, and orbits of a test particle 
moving under the pseudo-Newtonian potential, respectively. 
The run time is noted at the top-right corner in units of 



Q^ 1 , while the number of SPH particles are indicated at the 
bottom-right corner. The star is tidally disrupted before it 
passes through the first pericenter at t = 0.5. Afterward, the 
stellar debris expands along the orbit of the test particle as 
shown in panel (a) and (b). It is stretched towards the second 
pericenter as shown in panel (c). At the second pericenter, 
the stretched debris weakly intersects with the remnant of 
the debris in panel (d) and then reaches the second apocen- 
ter while continuing to expand, in panel (e). When it reaches 
the third pericenter in panel (f), two stretched orbits clearly 
cross over. The orbital energy is significantly dissipated by 
the shock from orbit crossing between the two stretched de- 
bris streams. This is consistent with the decrease of specific 
energy in panel (a) of Fig. 1111 The stellar debris is therefore 
rapidly circularized as shown in panels (g) and (h). Finally, 
in panel (i), a disk like structure is formed around the black 
hole sufficiently inside the tidal disruption radius. Note that 
the number of SPH particles are slightly reduced from panel 
(e) to panel (i) , since a very small fraction of the total num- 
ber of SPH particles is accreted onto the black hole. The 
trajectory of a SPH particle picked out of this simulation 
for ^ t < 4 is drawn in the solid line of Fig. [3] The size of 
the accretion disk is in rough agreement with r c irc ~ 2.5r p . 

Although the pseudo-Newtonian potential does not ac- 
curately model the precession of near-circular orbits, we 
note that by the time the pseudo-Newtonian potential be- 
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Figure 11. Evolution of the specific energy and specific angular momentum in Model 2a. In panel (a), the solid line and dashed line 
represent the specific energy of SPH particles and that of a test particle et p = —0.497, respectively. The small inside panel shows the 
evolution of the specific internal energy, ei nt . They are normalized by et = GM/rt- Panel (b) shows the specific angular momentum 
normalized by hp = 0.714, which is the specific angular momentum of the initial stellar orbit. 



comes significantly inaccurate, the debris streams have en- 
tered a regime of frequent orbit crossings, guaranteeing fur- 
ther rapid circularization. Furthermore, at a qualitative level 
the pseudo-Newtonian potential probably overestimates the 
circularization timescale, as it underestimates the true rate 
of relativistic precession for eccentric orbits. However, the 
pseudo-Newtonian potential still plays a crucial role in or- 
bital circularization processes, because the motion of a test 
particle under the Newtonian potential take place in a closed 
path and therefore causes no orbit crossing. In order to test 
this, we have performed the simulation for Model 2b, which 
has the same simulation parameters as Model 2a except that 



the Newtonian potential is adopted. Fig. [13] shows sequen- 
tial snapshots of the surface density of stellar debris for 
Model 2b. It has the same format as Fig. [T2j We note that 
there is no significant evidence for orbital crossings during 
the timescale of the simulation, since SPH particles orbit 
around the black hole on fixed eccentric orbits. This is also 
confirmed by the fact that the percentage change of specific 
orbital energy from t=0 to t—4 for Model 2b is less than 
0.4%. 

Fig. [14] shows the time-dependence of the number of 
SPH particles inside the tidal disruption radius iVacc and its 
first derivative (the mass capture rate). After the tidal dis- 
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Figure 12. A sequence of snapshots of the tidal disruption process in Model 2a. They are from panel (a) to panel (i) in chronological 
order. Each panel shows a surface density projected on x-y plane in five orders of magnitude in a logarithmic scale for 0.6 ^ t ^ 4, where 
t is in units of Q* 1 . The black hole is set at the origin. The run time is annotated at the top-right corner, while the number of SPH 
particles are indicated at the bottom-right corner. The dashed circle and dotted line indicate the tidal disruption radius and the orbits 
of a test particle moving under the pseudo-Newtonian potential given by equation (|10|) . respectively. 



ruption of the star, its orbit passes through the first apocen- 
ter, going completely outside of the tidal disruption radius. 
The first peak of Af aC c in panel (a) comes when the debris 
streams pass from the first apocenter to the second apoc- 
enter via the second pericenter, and the stretched debris 
re-enters the tidal disruption radius. Part of it exits once 
more, but the fractional remaining part is still inside the 
tidal disruption radius. A sequence of these events can be 
seen in panels (c)-(e) of Fig. [T2j The second peak of N acc 
forms in panel (b) . The stretched debris returns again to the 



tidal disruption radius, moving toward the third pericenter. 
Afterwards, most of debris circularizes and remains inside 
of the tidal disruption radius. 

Panel (b) of Fig. [14] shows the rate of mass being cap- 
tured inside the tidal disruption radius. The first three peaks 
are formed as stellar debris passes in and out the tidal dis- 
ruption radius, while the final peak shows the mass transfer 
rate to the accretion disk around the black hole. The mass 
transfer rate is also clipped from panel (b) to Fig. \TE\ as the 
mass accretion rate normalized by 1 Mq yr _1 . It is clear 
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Figure 13. Same format as Model 2a but for Model 2b. 



that the mass accretion rate has completely deviated from 
not only the canonical t _5//3 law but also the accretion rate 
variation of the other highly eccentric TDEs ( see Panel (b) 
of Figs. [9l and UOp . The resultant accretion rate is more than 
five orders of magnitude higher than the Eddington rate (see 
equation [17]). This shows that the accretion flow is extremely 
supercritical in the case of moderately eccentric TDEs. 

Since the viscous timescale measured at r c irc is esti- 
mated to be tvis ~ 3.5 x 10 3 (0.1/ass)^^ 1 where ass is the 
Shakura-Sunyaev viscosity parameter, the accretion flow is 
clearly in a non- steady state. Here, we assume that the ac- 
cretion disk is a geometrically thick : r c i rc /H ~ 1, where 
H is the disk scale height. Although the fate of the circu- 
larized debris is unclear because of much shorter simulation 



time than the viscous timescale, the super-Ed dington accre- 
tion flow will likely drive a powerful outflow (|Ohsuga et al.l 
120051 ) as it becomes radiation-pressure dominated. This may 
increas e the optical luminosity of t he flare by orders of mag- 
nitude (|Strubbe &; Quataerdl2009h . It is even possible that 
a radiation-pressure supported envelope could be formed 
(|Loeb fe Ulmerlll997l ). 



4 SOURCES OF ECCENTRIC TDES 

Although the canonical source for TDEs is two-body scat- 
tering at par sec scales, a variety of alternate mechanisms 
exist. In this section we briefly review these, and highlight 
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Figure 14. The number of SPH particles N acc entering inside the tidal disruption radius over ^ t ^ 4 in panel (a) and its corresponding 
mass capture rate, which is measured by the number of particles entering inside the tidal disruption radius per unit time. In panel (b), 
the mass capture rate is normalized by Mq/D" 1 , which corresponds to 8.9 x 10 5 MEdd for a 1O 6 M0 black hole with mass-to-energy 
conversion efficiency 77 = 0.1. 



those which could supply stars at eccentricities e < e cr it to 
their central SMBH. 

• Resonant relaxation: at small spatial scales close to a 
SMBH, the approximately Keplerian potential causes stel- 
lar angular momentum to relax in a correlated, coherent 
way (as opposed to the uncorrelated changes from two- 
body relaxation). It has been suggested that the rapid 
changes in orbital eccentricity produced by resonant re- 
laxation (RR) may en hance the TDE rate at small scales 
(jRauch fc Inga lls 1998); however, recent work has indicated 
that the complicated interplay between RR and general rel- 
ativistic precession gives rise to a "Schwarzschild barrier," 
which quenches RR for stars with small semi-major axis 
(|Merritt et al.l 1201 lh . RR is therefore unlikely to produce 
TDEs with e < e cr it. 

• Nuclear triaxiality: conversely, the triaxial, non- 
Keplerian potential of a surrounding star cluster results 
in fewer constants of motion for stellar orbits; in partic- 
ular, angular momentum is not conserved and some stars 
can wande r arbitrarily close to th e SMBH, eventually being 
disrupted (|Merritt &; Poonll2004l ). However, this effect only 
arises at large radii where the potential from the star cluster 
can induce a significant triaxial correction to the SMBH po- 
tential, and is therefore unlikely to produce low eccentricity 
TDEs. 

• A nuclear stellar disk: there is a rotating stellar disk 
comp osed of young mass ive stars at the center of the Milky 
Way (jBartko et alJl20ld ). These stars orbit with moderate 
eccentricity around the supermassive black hole Sgr A* . The 



nonresonant relaxation timescale drops as a star cluster be- 
comes flattened and disk-like, due to greater encounter rates 
between stars; simple estimates of the two-body eccentric- 
ity relaxation timescale in the Milky Way's stellar disk 
5000M©) find values t r i x - 4 x 10 8 yrs (|Kocsis k, Tremaind 
l201ll ). Assuming a typical stellar mass of 1OM gives a disk 
star TDE rate of ~ 10 -6 yr -1 , implying that if analogous 
disks are common in other galaxies, they could contribute 
nontrivially to the total TDE rate. However, in order to pro- 
duce TDEs with e < e cr it, disk stars would need to be scat- 
tered in from ~ mpc scales, implying that more compact 
disks than the Milky Way's are needed to produce signifi- 
cantly eccentric TDEs. 

• Binary SMBHs: A binary SMBH system will, for a 
period of ~ 10 5 yr to ~ 10 6 yr, see a TDE rate enhance- 
ment up to 10 _1 yr -1 from a combination of chaotic orbital 
evolu tion and Kozai cycles (|lvanov et al.l 120051 : 1 Chen et al.l 
l2009h . Recent work suggests that chaotic orbits in particu- 
lar are the dominant contribution to the rate enhancement 
(jChen et al.ll201ll : IWegg fc Nate Bodd lioilh . Most of the 
stellar flux originates from within one order of magnitude 
of the semimajor axis of the binary; since the inspiral of bi- 
nary SMBHs stalls when the lower-mass secondary reaches 
a radius inside of which is contained its own mass in stars, a 
low binary mass ratio q will pr oduce a flux of les s eccentric 
TDEs. In particular, Fig. 17 of IChen et all (|201lh indicates 
that for q = 1/81 and a primary black hole mass of 10 7 M©, 
some stellar flux into the primary's loss cone originates from 
spatial scales ^ 10 15 cm, which implies e < e cr it even for 
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P = 1 events. However, the chaotic orbits which dominate 
the r ate enhancement produ ce TDEs sampling a wide range 
of /3 (|Merritt &; Poon| [2004). so the situation is even more 
favorable. Binary SMBHs on a eccentric orbit appear to pro- 
duce even more TDEs from small apocenters. The "harden- 
ing radius" (the orbital radius within which the stellar cusp 
of the primary, larger SMBH contains the mass in stars of 
the secondary, smaller SMBH) effectively sets the scale at 
which stellar disruption rates are enhanced; to produce ec- 
centric TDEs requires relatively low mass ratios (< 1/50) 
which could in some cases be characterized as SMBH-IMBH 
inspirals. Because binary SMB Hs are expected to provide up 
to 10% of the total TDE rate (jWegg fc Nate Boddlioilh . it 
is important to note that a subset of these events may devi- 
ate from canonical "parabolic" light curves. 

• SMBH recoil: the gravitational wave recoil accompa- 
nying a SMBH merger will scramble the orbits of surround- 
ing stars and partially refill the kicked merger remnant's 
loss cone, briefly increasing the TDE rate. If one considers 
the bulk of the stellar population surrounding the merging 
binary black holes, the excavation of a "binary loss cone" 
results in too few stars at small separation for the result- 
i ng burst of TDEs to involve any events with e < e cr it 
(jStone fe Loebll201lh . However, at low binary mass ratios, it 
is possible for mean motion resonances to mig rate stars to 
small semimajor axes dur i ng the binary inspiral (|Schnittmanl 
l20ld : ISeto fc MutolbOloL l201lh : any TDEs produced from 
this stellar subpopulation would have e < e cr it, and would 
be described by this paper. 

• Binary separation: Recent theoretical studies on 
rates of tidal separation of binary stars by SMBHs suggest 
that a significant fraction of tidal disruption flares may occur 
from stars approaching the black hole from subpa rsec scales 
(jAmaro-Seoane et al.l l2Qlll : iBromlev et~all l2012h . When a 
binary star passes sufficiently close to the black hole to be 
tidally separated (without immediate disruption of either 
component), one star becomes tightly bound to the black 
hole while the other is flung away at a high speed. The sub- 
sequent orbital evolution of the bound star will represent 
a competition between gravitational radiation and stellar 
relaxation processes. If gravitational radiation dominates, 
the orbit will circularize and spiral in, likely leading to a 
phase of steady mass transfer that is unlike the eccentric 
and violent disruptions which we have simulated in this 
paper. On the other hand, if stellar relaxational processes 
dominate, the star can eventually diffuse into the loss cone 
with a nonzero eccentricity that is still significantly smaller 
than that expected for TDEs generated by two-body relax- 
ation at parsec scales. For a bound star with r p > r t and 
semimajor axis a, living in a stellar cusp with density pro- 
file p(r) oc r~ b , the eccentricity dividing gravitational-wave- 
dominated evolution from r elaxation-dominated evolution is 
(jAmaro-Seoane et aHl2Qllh 

e gw « 1 - 0.016 x (8 x 10- 4 ) (2b - 3)/5 (3 - b) 2/5 (24) 



( M BH ^ 

V lo7M © J 



(8-b)/5 



(lOOO Au) 



(2b-ll)/5 



In practice, requiring e gw > e cr it for a < r t /(l — e cr it) is 
only possible for steep stellar cusps b ~ 2 and very low 
mass SMBHs (Mbh < 10 5-5 M©). A further complication is 
the presence of the Schwarzschild barrier, which may pre- 
vent stars with long gravitational wave inspiral timescales 



from entering the loss cone. Stars bound to a SMBH from 
binary separation events can become TDEs with lower ec- 
centricity than in the canonical scenario, but equation (|24[) 
implies their eccentricity will still tend to be greater than 
e cr it- On the other hand, equation ([24]) was derived assum- 
ing the Fokker-Planck diffusion limit for two-body relax- 
ation, and anomalou s diffusion from strong scattering events 
(|Bar-Or et al.l l2012h may allow violent tidal disruption of 
separated binaries even from eccentricities e < e gw . 

Of the possibilities we have considered, two stand out as 
particularly promising: TDEs generated as binary SMBHs 
harden and stall, and TDEs generated during the coales- 
cence of binary SMBHs with a low mass ratio by stars that 
were brought inward via mean motion resonances. The first 
possibility may account for a significant fraction of the total 
TDE rate; the second possibility will be much less intrinsi- 
cally common, but could serve as an electromagnetic coun- 
terpart to an eLISA/NGCEl-band gravitational wave signal. 
It is also possible that strong scatterings could lead to ec- 
centric TDEs from the bound stars produced by tidal sep- 
aration of stellar binaries. Other possible enhancements to 
theoretical TDE rates, such as perturbations from infalling 
giant molecular clouds, likely occur at too large of a spatial 
scale to create eccentric TDEs. 



5 SUMMARY & DISCUSSION 

We have carried out numerical simulations of mass fallback 
and accretion processes around a SMBH. Specifically, we 
have examined the tidal disruption of a star on a bound or- 
bit, considering relativistic effects with a pseudo-Newtonian 
potential. Using a polytropic gas sphere as our initial con- 
ditions, we have considered both parabolic orbits (e = 1.0) 
and eccentric ones (e = 0.8 and 0.98), varying the penetra- 
tion factor f3 as well. 

We have found that a non-steady, non-axisymmetric ac- 
cretion disk is formed around the black hole in the case of 
e=0.8 and /3 = 5. The formation of an accretion disk oc- 
curs as follows: a segment of the stellar debris returning 
to pericenter, and a different one exiting pericenter, inter- 
sect each other due to relativistic precession. The orbital 
energy is then dissipated by shocks due to orbit crossing. 
Since the orbital angular momentum of the stellar debris is 
conserved before and after the tidal disruption, the debris 
orbits are rapidly circularized during a few orbit crossings. 
This shows that the initial size of the accretion disk is only 
determined by the orbital angular momentum of the initial 
stellar orbit. In our simulations, the circularization radius 
is estimated aS Tcirc ~ 2.5r p where r p is the pericenter dis- 
tance. Furthermore, the expected accretion rate is extremely 
super-Eddington. 

The striking difference between moderate-eccentricity 
simulations with and without relativistic precession high- 
lights the importance of general relativistic effects for de- 
bris circularization. Specifically, very little energy dissipa- 
tion was seen in Model 2b (Newtonian) of Table 1, while its 
pseudo-Newtonian equivalent, Model 2a, saw rapid accre- 
tion disk formation. For reasons of computational cost, we 
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did not follow the longer timescales required for debris cir- 
cularization in other models, but in future work we hope 
to investigate if these results generalize to higher eccen- 
tricity and parabolic TDEs. Sufficiently rapid SMBH spin 
could delay disk formation, as Lense-Thirring torques will 
precess the orbital planes of individual debris streams and 
limit or prevent orbit-crossings. If a disk is able to eventu- 
ally form, its luminosity will b e periodically modula ted by 
Lense-Thirring precession fe.g. JStone fc Loebll2012ah . Even 
considering this complication, the energy dissipation process 
is crucial for determining the formation and structure of an 
accretion disk. In these exploratory simulations, we have 
adopted the polytropic equation of state instead of solving 
an energy equation. In a subsequent paper, we will study the 
detailed disk formation and structure by solving more real- 
istic energy equations with and without a radiative cooling 
term. 

Eccentric TDEs are a subpopulation of all TDEs, but a 
potentially interesting one. The two most promising means 
of producing them are in the dynamical friction stage of bi- 
nary SMBHs, and immediately following gravitational wave- 
driven black hole coalescence. Both of these are especially 
interesting subsets of TDEs: the former occurs during a time 
of greatly enhanced TDE rates, when it may be possible for 
a single galaxy to produce multiple TDEs in ~ 10 yr, while 
the latter would serve as a delayed electromagnetic coun- 
terpart to a low-frequency gravitational wave signal. It is 
also possible that tidal separation of binary stars could pro- 
duce TDEs of interestingly low-to-moderate eccentricity, if a 
large population of Mbh ^ 1O 5 ' 5 M SMBHs reside in steep 
stellar cusps. Eccentric TDEs are also of interest because of 
the relatively short delay time between disruption and disk 
formation, allowing numerical simulations to bypass costly 
apocenter passages and directly approach open problems in 
debris cir culariz a tion. T he pseudo-Newtonian potential pro- 
posed by IWegd {2012) is, however, not applicable for the 
very low eccentricities of fully circularized gas because it 
was derived with the assumption of small binding energy. 
Therefore, we hope to employ post-Newtonian corrections 
in a subsequent paper. 

For eccentric stellar orbits, there is a critical value of the 
orbital eccentricity below which all stellar debris is bounded 
to the black hole. It can be seen from our simulations that 
the critical eccentricity is slightly lower than expected from 
our analytical predictions. This might be because of the 
effects of stellar structure on the tidal disruption. There 
are three important implications for lightcurves of eccentric 
TDEs: 

• The mass fallback rate will not asymptote to M oc £ _5//3 
at late times but be finite with cut-off time t = t max - 

• The lack of unbound debris will eliminate obser- 
vational signatures assoc iated with emission lines (e.g., 
IStrubbe fe Quataertll2009h . 

• A larger amount of mass will return to pericenter in a 
much shorter time than in the standard parabolic picture, 
considerably increasing the ratio of M/M^dd for the subse- 
quent flare. 

Even for eccentricities e > e cr it , the center of the differential 
mass distribution will shift in a negative direction, provid- 
ing weaker versions of the above effects. These signatures 
should be searched for when large samples of TDE candi- 
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Figure 15. Mass accretion rate normalized by 1 solar mass per 
year, which corresponds to the final peak in panel (b) of Fig. 1141 

dates from next generation optical transient surveys such 
as the Large Synoptic Survey Telescope^ and next genera- 
tion all-sky X-ray surveys such as extended Roentgen Survey 
with an Imaging Telescope Arrajjf] become available. 
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